#######################################################################################################
## Replication Script                                                                                ##    
## State Capacity, Refugee Reception, and Attitudes towards Refugees                                 ## 
## Study 1: Attitudes towards Ukrainians Refugees in CEE                                             ##  
#######################################################################################################


#######################################################################################################
#### packages ####
library(rio)
library(dplyr)
library(tidyr)
library(stringr)
library(ggplot2)
library(ggtext)
library(ggpubr)
library(patchwork)
library(modelsummary)
library(kableExtra)
library(polycor)
library(sandwich)
library(clubSandwich)
library(lmtest)
library(car)
library(Amelia)
library(mice)
library(howManyImputations)

set.seed(0987)

#######################################################################################################

#######################################################################################################
#### sessionInfo() ####

# R version 4.4.0 (2024-04-24)
# Platform: aarch64-apple-darwin20
# Running under: macOS Sonoma 14.5

# Matrix products: default
# BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib  
# LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.0

#locale:
#  [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

# time zone: America/Chicago
# tzcode source: internal

# attached base packages:
#  [1] stats     graphics  grDevices utils     datasets  methods   base   

# other attached packages:
#  [1] howManyImputations_0.2.5 mice_3.16.0              Amelia_1.8.2             Rcpp_1.0.13          
#  [5] car_3.1-2                carData_3.0-5            lmtest_0.9-40            zoo_1.8-14  
#  [9] clubSandwich_0.5.11      sandwich_3.1-1           polycor_0.8-1            kableExtra_1.4.0      
# [13] modelsummary_2.1.1       patchwork_1.2.0          ggpubr_0.6.0             ggtext_0.1.2  
# [17] ggplot2_3.5.2            stringr_1.5.1            tidyr_1.3.1              dplyr_1.1.4  
# [21] rio_1.2.0                

# loaded via a namespace (and not attached):
# [1] rlang_1.1.6         magrittr_2.0.3      compiler_4.4.0      systemfonts_1.1.0   vctrs_0.6.5         pkgconfig_2.0.3   
# [7] shape_1.4.6.1       crayon_1.5.3        fastmap_1.2.0       backports_1.5.0     labeling_0.4.3      effectsize_0.8.9  
# [13] rmarkdown_2.27      markdown_1.13       tzdb_0.4.0          haven_2.5.4         nloptr_2.1.1        purrr_1.0.2      
#[19] xfun_0.46           glmnet_4.1-8        jomo_2.7-6          tinytable_0.3.0     pan_1.9             broom_1.0.6       
#[25] parallel_4.4.0      R6_2.6.1            tables_0.9.28       stringi_1.8.4       RColorBrewer_1.1-3  parallelly_1.37.1 
#[31] boot_1.3-30         rpart_4.1.23        estimability_1.5.1  iterators_1.0.14    knitr_1.48          future.apply_1.11.2
#[37] parameters_0.24.1   R.utils_2.12.3      readr_2.1.5         Matrix_1.7-0        splines_4.4.0       nnet_7.3-19       
#[43] tidyselect_1.2.1    rstudioapi_0.16.0   abind_1.4-5         codetools_0.2-20    listenv_0.9.1       admisc_0.38       
#[49] lattice_0.22-6      tibble_3.3.0        bayestestR_0.15.2   withr_3.0.2         coda_0.19-4.1       evaluate_0.24.0   
#[55] foreign_0.8-87      future_1.33.2       survival_3.7-0      xml2_1.3.6          pillar_1.11.0       checkmate_2.3.1   
#[61] foreach_1.5.2       insight_1.0.2       generics_0.1.3      hms_1.1.3           commonmark_1.9.1    scales_1.4.0      
#[67] minqa_1.2.7         xtable_1.8-4        globals_0.16.3      glue_1.8.0          emmeans_1.10.3      tools_4.4.0       
#[73] data.table_1.16.0   lme4_1.1-35.5       ggsignif_0.6.4      forcats_1.0.0       mvtnorm_1.2-5       grid_4.4.0        
#[79] datawizard_1.0.0    nlme_3.1-165        performance_0.12.2  cli_3.6.5           fansi_1.0.6         viridisLite_0.4.2 
#[85] svglite_2.1.3       gtable_0.3.6        R.methodsS3_1.8.2   rstatix_0.7.2       digest_0.6.36       farver_2.1.2      
#[91] htmltools_0.5.8.1   R.oo_1.26.0         lifecycle_1.0.4     mitml_0.4-5         gridtext_0.1.5      MASS_7.3-61  
#######################################################################################################

#######################################################################################################
#### Figure 1 Data Import ####

eb1 <- rio::import("EB_97_5.dta") 
eb2 <- rio::import("EB_98_2.dta")
eb3 <- rio::import("EB_99_4.dta")
eb4 <- rio::import("EB_100_2.dta")
eb5 <- rio::import("EB_101_3.dta")

#######################################################################################################

#######################################################################################################
#### Figure 1 Data Cleaning ####

eb1 <- eb1[eb1$isocntry=="CZ" | eb1$isocntry=="HU" | eb1$isocntry=="PL" | eb1$isocntry=="SK", ]
eb2 <- eb2[eb2$isocntry=="CZ" | eb2$isocntry=="HU" | eb2$isocntry=="PL" | eb2$isocntry=="SK", ]
eb3 <- eb3[eb3$isocntry=="CZ" | eb3$isocntry=="HU" | eb3$isocntry=="PL" | eb3$isocntry=="SK", ]
eb4 <- eb4[eb4$isocntry=="CZ" | eb4$isocntry=="HU" | eb4$isocntry=="PL" | eb4$isocntry=="SK", ]
eb5 <- eb5[eb5$isocntry=="CZ" | eb5$isocntry=="HU" | eb5$isocntry=="PL" | eb5$isocntry=="SK", ]

eb1 <- select(eb1, uniqid, isocntry, accept_ukr_ref = qe2_5)
eb2 <- select(eb2, uniqid, isocntry, accept_ukr_ref = qe2_5)
eb3 <- select(eb3, uniqid, isocntry, accept_ukr_ref = qe2_5)
eb4 <- select(eb4, uniqid, isocntry, accept_ukr_ref = qd2_5)
eb5 <- select(eb5, uniqid, isocntry, accept_ukr_ref = qd2_5)

eb1$year <- "July 2022"
eb2$year <- "Feb. 2023"
eb3$year <- "June 2023"
eb4$year <- "Nov. 2023"
eb5$year <- "May 2024"

eb <- rbind(eb1,eb2,eb3,eb4,eb5)

table(eb$accept_ukr_ref)

eb <- eb %>%
  mutate(accept_ukr_ref = case_when(
    accept_ukr_ref == 1 ~ 4,
    accept_ukr_ref == 2 ~ 3,
    accept_ukr_ref == 3 ~ 2,
    accept_ukr_ref == 4 ~ 1,
    accept_ukr_ref == 5 ~ NA_real_,
    TRUE ~ accept_ukr_ref
  ))

table(eb$accept_ukr_ref)

eb_agg <- eb %>%
  group_by(isocntry, year) %>%
  summarise(across(starts_with("accept_ukr_ref"), list(Mean = ~mean(., na.rm = TRUE), CI_Lower = ~t.test(., na.rm = TRUE)$conf.int[1], CI_Upper = ~t.test(., na.rm = TRUE)$conf.int[2])))

eb_percent <- eb %>%
  group_by(isocntry, year) %>%
  summarise(
    n_total = n(),
    n_low = sum(accept_ukr_ref %in% 3:4, na.rm = TRUE),
    pct_low = 100 * n_low / n_total
  ) %>%
  ungroup()

eb_percent <- eb_percent %>%
  mutate(year = factor(year, levels = c(
    "July 2022",
    "Feb. 2023",
    "June 2023",
    "Nov. 2023",
    "May 2024"
  )))

eb_percent$isocntry[eb_percent$isocntry == "CZ"] <- "Czech Rep."
eb_percent$isocntry[eb_percent$isocntry == "HU"] <- "Hungary"
eb_percent$isocntry[eb_percent$isocntry == "PL"] <- "Poland"
eb_percent$isocntry[eb_percent$isocntry == "SK"] <- "Slovakia"

#######################################################################################################

#######################################################################################################
#### Figure 1 Plot ####

eb_attitude_plot <- ggplot(eb_percent, aes(x = year, y = pct_low, color = isocntry, group = isocntry)) +
  geom_line(size = 0.6) +
  geom_point(size = 1.25) +
  ylim(50,100) +
  scale_color_manual(
    values = c(
      "Czech Rep." = "#999999",
      "Slovakia" = "#E69F00",
      "Poland" = "#56B4E9",  
      "Hungary" = "#009E73")) +
  labs(
    title = "",
    x = "Survey Date",
    y = "% Agree",
    color = "Country"
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    strip.background = element_rect(fill = "#f0f0f0", color = "gray70"),
    strip.text = element_text(size = 10),
    legend.position = "bottom", legend.title = element_blank())

eb_attitude_plot

# ggsave("Fig1_EB_Attitude_Plot.jpeg", plot = eb_attitude_plot, width=4 ,height = 3.5, dpi=600)

rm(list = ls())

#######################################################################################################

#######################################################################################################
#### Study 1 Data Import + Data Cleaning ####

data <- rio::import("Hlatky_UKR_Ref_Study1_Survey_Data.csv")

#reverse code grievance/concerns so that higher values indicate stronger grievance ### 


data <- data %>%
  mutate(
    Q16A  = case_when(Q16A  == 1 ~ 4,
                      Q16A  == 2 ~ 3,
                      Q16A  == 3 ~ 2,
                      Q16A  == 4 ~ 1,
                      TRUE       ~ Q16A),
    
    Q16C = case_when(Q16C  == 1 ~ 4,
                     Q16C  == 2 ~ 3,
                     Q16C  == 3 ~ 2,
                     Q16C  == 4 ~ 1,
                     TRUE       ~ Q16C),
    
    Q16D = case_when(Q16D  == 1 ~ 4,
                     Q16D  == 2 ~ 3,
                     Q16D  == 3 ~ 2,
                     Q16D  == 4 ~ 1,
                     TRUE       ~ Q16D),
    
    Q16E = case_when(Q16E  == 1 ~ 4,
                     Q16E  == 2 ~ 3,
                     Q16E  == 3 ~ 2,
                     Q16E  == 4 ~ 1,
                     TRUE       ~ Q16E),
    
    Q16H  = case_when(Q16H  == 1 ~ 4,
                      Q16H  == 2 ~ 3,
                      Q16H  == 3 ~ 2,
                      Q16H  == 4 ~ 1,
                      TRUE       ~ Q16H),
    
    Q16J  = case_when(Q16J  == 1 ~ 4,
                      Q16J  == 2 ~ 3,
                      Q16J  == 3 ~ 2,
                      Q16J  == 4 ~ 1,
                      TRUE       ~ Q16J)
  )
#######################################################################################################

#######################################################################################################
#### Study 1: Table A1: Study 1 Descriptive Statistics  ####

datasummary(Q2 + Q16A + Q16B + Q16C + Q16D + Q16E + Q16G + Q16H + Q16I + Q16J + Q21 + as.factor(Q1) + as.factor(Q3) + as.factor(Q4a) + as.factor(Q7_3) + as.factor(Q7_5) + as.factor(Q11_1) + as.factor(Q11_2) + as.factor(Q11_3) + as.factor(Q11_4) + as.factor(Q11_7) + as.factor(Q12_1) + as.factor(Q12_2) + as.factor(Q12_3) + as.factor(Q12_4) + as.factor(Q14) + as.factor(Q8) + as.factor(COUNTRY) ~ Mean + SD + Min + Max + N, data=data, dinm=FALSE, output = "latex")

#######################################################################################################

#######################################################################################################
#### Study 1: Figure 2: Concerns with Ukrainian Refugees in Central Europe  ####

data_long <- data %>%
  pivot_longer(cols = c(Q16A, Q16H, Q16G, Q16J, Q16I), 
               names_to = "item", values_to = "response") %>%
  mutate(
    grievance = case_when(
      item %in% c("Q16A", "Q16H", "Q16G", "Q16J", "Q16I") ~ 
        ifelse(is.na(response), 0, ifelse(response >= 3, 1, 0)),
      TRUE ~ 0
    )
  )

data_long$item <- factor(data_long$item, levels = c("Q16H", "Q16I", "Q16A", "Q16G", "Q16J"))

summary_data <- data_long %>%
  group_by(country, item) %>%
  summarise(
    n_total = n(),                                
    n_grievance = sum(grievance == 1),
    grievance_percent = (n_grievance / n_total) * 100,
    se = sqrt((grievance_percent / 100) * (1 - grievance_percent / 100) / n_total),
    lower = grievance_percent - 1.96 * se * 100,
    upper = grievance_percent + 1.96 * se * 100,
    .groups = "drop"
  )

labels <- c(
  Q16A = "Reduce Support",
  Q16H = "Labor Market",
  Q16G = "Cultural",
  Q16J = "Safety",
  Q16I = "Economy"
)

grievance_plot <- ggplot(summary_data, aes(x = country, y = grievance_percent, color=country)) +
  geom_point(size = 1.25) +
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0) +
  scale_color_manual(values = c(
    "CZ" = "#999999",
    "SK" = "#E69F00",
    "PL" = "#56B4E9",  
    "HU" = "#009E73"
  ), guide = "none") +
  facet_wrap(~ item, labeller = as_labeller(labels), ncol = 5) +
  labs(
    y = "% Expressing Concern",
    x = "Country",
  ) +
  theme_minimal() +
  theme(
    strip.background = element_rect(fill = "#f0f0f0", color = "gray70"),
    strip.text = element_text(size = 10),
    panel.border = element_rect(color = "black", fill = NA, linewidth = 0.25)
  ) + ylim(0,80)

grievance_plot

# ggsave("Fig2_Concerns_UKR_Ref.jpeg", plot = grievance_plot, width=6.75 ,height = 4.5, dpi=600)

#######################################################################################################

#######################################################################################################
#### Study 1: Figure A1: Concerns with Ukrainian Refugees in Central Europe (Original Scale) ####


data_long <- data %>%
  pivot_longer(cols = c(Q16A, Q16H, Q16G, Q16J, Q16I), 
               names_to = "item", values_to = "response")

data_long$item <- factor(data_long$item, levels = c("Q16H", "Q16I", "Q16A", "Q16G", "Q16J"))

summary_data <- data_long %>%
  filter(!is.na(response)) %>%  # Remove NAs before calculating statistics
  group_by(country, item) %>%
  summarise(
    n_total = n(),                                
    mean_response = mean(response),
    se = sd(response) / sqrt(n_total),
    lower = mean_response - 1.96 * se,
    upper = mean_response + 1.96 * se,
    .groups = "drop"
  )



labels <- c(
  Q16A = "Reduce Support",
  Q16H = "Labor Market",
  Q16G = "Cultural",
  Q16J = "Safety",
  Q16I = "Economy"
)

grievance_plot_alt <- ggplot(summary_data, aes(x = country, y = mean_response, color = country)) +
  geom_point(size = 1.25) +
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0) +
  scale_color_manual(values = c(
    "CZ" = "#999999",
    "SK" = "#E69F00",
    "PL" = "#56B4E9",  
    "HU" = "#009E73"
  ), guide = "none") +
  facet_wrap(~ item, labeller = as_labeller(labels), ncol = 5) +
  labs(
    y = "Mean Response",
    x = "Country"
  ) +
  theme_minimal() +
  theme(
    strip.background = element_rect(fill = "#f0f0f0", color = "gray70"),
    strip.text = element_text(size = 10),
    panel.border = element_rect(color = "black", fill = NA, linewidth = 0.25)
  )

grievance_plot_alt

# ggsave("FigA1_Concerns_UKR_Ref_alt.jpeg", plot = grievance_plot_alt, width=6.75 ,height = 4.5, dpi=600)

#######################################################################################################

#######################################################################################################
#### Study 1: Table A2: Correlations of Predictors and Outcomes ####

items <- c("Q11_7", "Q16A", "Q16H", "Q16G", "Q16J","Q16I")

data_cor <- data[, items]

cor_matrix <- hetcor(data_cor)  
print(cor_matrix)                 

data_cor <- data_cor %>%
  mutate(
    Q11_7 = ordered(Q11_7, levels = 0:1),
    across(c(Q16A, Q16H, Q16G, Q16J, Q16I), ~ordered(., levels = 1:4))
  )

cor_matrix <- hetcor(data_cor)  

cor_df <- cor_matrix$correlations
var_order <- c("Q11_7", "Q16H", "Q16I", "Q16A", "Q16G", "Q16J")
cor_df <- round(cor_matrix$correlations, 4)
cor_df <- cor_df[var_order, var_order]

cor_df[upper.tri(cor_df)] <- ""

print(cor_df)

kable(cor_df, "latex",
      caption = "Correlations of Predictors and Outcomes",
      booktabs = TRUE, digits = 4, na = "") %>%
  kable_styling(latex_options = c("hold_position"), full_width = FALSE) %>%
  row_spec(0, bold = TRUE) %>%
  column_spec(1, bold = TRUE) %>%
  footnote(general = "Polychoric correlations.")

#######################################################################################################

#######################################################################################################
#### Study 1: Table 1, Table A3: Refugee-Related Concerns and Preferences over Ukrainian Refugee Admissions ####

model_allow_1 <- glm(Q11_7 ~ Q16A + Q16G + Q16H + Q16I + Q16J, data=data, family = "binomial")
summary(model_allow_1)

model_allow_2 <- glm(Q11_7 ~ Q16A + Q16G + Q16H + Q16I + Q16J + as.factor(Q4a) + as.factor(Q7_3) + as.factor(Q7_5) + as.factor(Q14)+ as.factor(Q1)+Q2+ as.factor(Q3) + Q21 +as.factor(COUNTRY) + as.factor(Q8), data=data, family="binomial")
summary(model_allow_2)

model_allow_3 <- glm(Q11_7 ~ Q16A + Q16G + Q16H + Q16I + Q16J + as.factor(Q4a) + as.factor(Q7_3)  + as.factor(Q7_5) + as.factor(Q14)+ as.factor(Q1)+Q2+ as.factor(Q3) + Q21 + as.factor(Q8), data=data[data$country=="CZ", ], family="binomial")
summary(model_allow_3)

model_allow_4 <- glm(Q11_7 ~ Q16A + Q16G + Q16H + Q16I + Q16J + as.factor(Q4a) + as.factor(Q7_3) + as.factor(Q7_5) + as.factor(Q14)+ as.factor(Q1)+Q2+ as.factor(Q3) + Q21 + as.factor(Q8), data=data[data$country=="HU", ], family="binomial")
summary(model_allow_4)

model_allow_5 <- glm(Q11_7 ~ Q16A + Q16G + Q16H + Q16I + Q16J + as.factor(Q4a)+ as.factor(Q7_3) + as.factor(Q7_5) + as.factor(Q14)+ as.factor(Q1)+Q2+ as.factor(Q3) + Q21 + as.factor(Q8), data=data[data$country=="PL", ], family="binomial")
summary(model_allow_5)

model_allow_6 <- glm(Q11_7 ~ Q16A + Q16G + Q16H + Q16I + Q16J + as.factor(Q4a) +  as.factor(Q7_3) + as.factor(Q7_5) + as.factor(Q14)+ as.factor(Q1)+Q2+ as.factor(Q3) + Q21 + as.factor(Q8), data=data[data$country=="SK", ], family="binomial")
summary(model_allow_6)


table1_models <- list(model_allow_1,model_allow_2,model_allow_3,model_allow_4,model_allow_5,model_allow_6)

vcov_list_table1 <- list(
  vcovCR(model_allow_1, cluster = data$COUNTRY, type = "CR2"),
  vcovCR(model_allow_2, cluster = data$COUNTRY, type = "CR2"),
  vcovHC(model_allow_3, type = "HC2"),
  vcovHC(model_allow_4, type = "HC2"),
  vcovHC(model_allow_5, type = "HC2"),
  vcovHC(model_allow_6, type = "HC2")
)

coef_map_full <- c("Q16H"="Labor Market", "Q16I"="Economy", "Q16A"="Reduce Support" , "Q16G"="Cultural Threat", "Q16J"="Safety", "as.factor(Q7_3)3"="Econ. Situation", "as.factor(Q7_5)5"="Govt. Eval.", "as.factor(Q4a)1"="Urban", "as.factor(Q14)2"="No Inappropriate Beh.", "as.factor(Q1)2"="Female", "Q2"="Age", "as.factor(Q3)2"= "Secondary, no leaving exam", "as.factor(Q3)3"= "Secondary, leaving exam", "as.factor(Q3)4"= "University","Q21"= "Household Financial Status", "as.factor(Q8)2" = "War Responsibility: Ukraine","as.factor(Q8)3" = "War Responsibility: US/NATO", "as.factor(COUNTRY)2"="CZ", "as.factor(COUNTRY)3"="PL", "as.factor(COUNTRY)4"="HU")


modelsummary(table1_models, output = "latex", coef_map = coef_map_full,  gof_omit = "AIC|RMSE", vcov = vcov_list_table1, stars = c('*' = .05,'**'=0.01))

#######################################################################################################
#### Study 1: Table A4: Variance Inflation Facotrs (VIF) for Table A3 Models ####

vif_model_1 <- vif(model_allow_1)
vif_model_2 <- vif(model_allow_2)
vif_model_3 <- vif(model_allow_3)
vif_model_4 <- vif(model_allow_4)
vif_model_5 <- vif(model_allow_5)
vif_model_6 <- vif(model_allow_6)

vif_model_1 <- as.matrix(vif_model_1) #since dif. number of dimensions


vif_df1a <- data.frame(
  Predictor = rownames(vif_model_1), 
  VIF = as.numeric(vif_model_1[, 1]), 
  model = "M1")

vif_df2a <- data.frame(
  Predictor = rownames(vif_model_2),
  VIF = as.numeric(vif_model_2[, 1]),
  model = "M2"
)
vif_df3a <- data.frame(
  Predictor = rownames(vif_model_3),
  VIF = as.numeric(vif_model_3[, 1]),
  model = "M3"
)
vif_df4a <- data.frame(
  Predictor = rownames(vif_model_4),
  VIF = as.numeric(vif_model_4[, 1]),
  model = "M4"
)
vif_df5a <- data.frame(
  Predictor = rownames(vif_model_5),
  VIF = as.numeric(vif_model_5[, 1]),
  model = "M5"
)
vif_df6a <- data.frame(
  Predictor = rownames(vif_model_6),
  VIF = as.numeric(vif_model_6[, 1]),
  model = "M6"
)


vif_all <- bind_rows(
  vif_df1a, vif_df2a, vif_df3a, vif_df4a, vif_df5a, vif_df6a)

vif_subset <- vif_all %>%
  filter(Predictor %in% c("Q16H", "Q16I", "Q16A", "Q16G", "Q16J")) %>%
  mutate(Predictor = factor(Predictor, levels = c("Q16H", "Q16I", "Q16A", "Q16G", "Q16J"))) %>%
  arrange(Predictor)

vif_wide <- vif_subset %>%
  pivot_wider(names_from = Predictor, values_from = VIF) %>%
  arrange(model)

vif_wide %>%
  kable(format = "latex", booktabs = TRUE, digits = 3,
        caption = "Variance Inflation Factors (VIF) for Table AX, Models") %>%
  kable_styling(latex_options = c("hold_position", "scale_down"))

rm(list = ls())

#######################################################################################################
#### Study 1: Table A5: Refugee-Related Concerns and Preferences over Ukrainian Refugee Admission (Multiple Imputation) ####

### re-import data + multiple imputation steps ### 
# data <- rio::import("Hlatky_UKR_Ref_Study1_Survey_Data.csv")

# reverse code grievance/concerns so that higher values indicate stronger 
# data <- data %>%
#  mutate(
#    Q16A  = case_when(Q16A  == 1 ~ 4,
#                      Q16A  == 2 ~ 3,
#                      Q16A  == 3 ~ 2,
#                      Q16A  == 4 ~ 1,
#                      TRUE       ~ Q16A),
#    
#    Q16C = case_when(Q16C  == 1 ~ 4,
#                     Q16C  == 2 ~ 3,
#                     Q16C  == 3 ~ 2,
#                     Q16C  == 4 ~ 1,
#                     TRUE       ~ Q16C),
#    
#    Q16D = case_when(Q16D  == 1 ~ 4,
#                     Q16D  == 2 ~ 3,
#                     Q16D  == 3 ~ 2,
#                     Q16D  == 4 ~ 1,
#                     TRUE       ~ Q16D),
#    
#    Q16E = case_when(Q16E  == 1 ~ 4,
#                     Q16E  == 2 ~ 3,
#                     Q16E  == 3 ~ 2,
#                     Q16E  == 4 ~ 1,
#                     TRUE       ~ Q16E),
#    
#    Q16H  = case_when(Q16H  == 1 ~ 4,
#                      Q16H  == 2 ~ 3,
#                      Q16H  == 3 ~ 2,
#                      Q16H  == 4 ~ 1,
#                      TRUE       ~ Q16H),
#    
#    Q16J  = case_when(Q16J  == 1 ~ 4,
#                      Q16J  == 2 ~ 3,
#                      Q16J  == 3 ~ 2,
#                      Q16J  == 4 ~ 1,
#                      TRUE       ~ Q16J)
#  )

# imputed_data <- amelia(data, idvars=c("qn","COUNTRY","country"), m=50, ords=c("Q11_7","Q11_1", "Q11_2", "Q11_3", "Q11_4", "Q16A", "Q16H", "Q16G", "Q16J","Q4a","Q14","Q1","Q3","Q8"))$imputations

# save(imputed_data, file = "Hlatky_UKR_Ref_Study1_Imputations.RData")

imputed_data <- rio::import("Hlatky_UKR_Ref_Study1_Imputations.RData")

mod <- list()

mod[['(1)']] <- lapply(imputed_data, function(x) glm(Q11_7 ~ Q16A + Q16G + Q16H + Q16I + Q16J, family = "binomial", x))
mod[['(1)']]  <- mice::pool(mod[['(1)']])



mod[['(2)']] <- lapply(imputed_data, function(x) glm(Q11_7 ~ Q16A + Q16G + Q16H + Q16I + Q16J + as.factor(Q4a) + as.factor(Q7_3) + as.factor(Q7_5) + as.factor(Q14)+ as.factor(Q1)+Q2+ as.factor(Q3) + Q21 +as.factor(COUNTRY) + as.factor(Q8), family = "binomial", x))
mod[['(2)']]  <- mice::pool(mod[['(2)']])
how_many_imputations(mod[['(2)']], cv = 0.05, alpha = 0.05)

mod[['(3)']] <- lapply(imputed_data, function(x) {
  glm(Q11_7 ~ Q16A + Q16G + Q16H + Q16I + Q16J + as.factor(Q4a) + as.factor(Q7_3)  + as.factor(Q7_5) + as.factor(Q14)+ as.factor(Q1)+Q2+ as.factor(Q3) + Q21 + as.factor(Q8),
      family = "binomial",
      data = subset(x, country == "CZ"),
      x = TRUE)
})
mod[['(3)']]  <- mice::pool(mod[['(3)']])
how_many_imputations(mod[['(3)']], cv = 0.05, alpha = 0.05)


mod[['(4)']] <- lapply(imputed_data, function(x) {
  glm(Q11_7 ~ Q16A + Q16G + Q16H + Q16I + Q16J + as.factor(Q4a) + as.factor(Q7_3) + as.factor(Q7_5) + as.factor(Q14)+ as.factor(Q1)+Q2+ as.factor(Q3) + Q21 + as.factor(Q8),
      family = "binomial",
      data = subset(x, country == "HU"),
      x = TRUE)
})
mod[['(4)']]  <- mice::pool(mod[['(4)']])
how_many_imputations(mod[['(4)']], cv = 0.05, alpha = 0.05)


mod[['(5)']] <- lapply(imputed_data, function(x) {
  glm(Q11_7 ~ Q16A + Q16G + Q16H + Q16I + Q16J + as.factor(Q4a) + as.factor(Q7_3) + as.factor(Q7_5) + as.factor(Q14)+ as.factor(Q1)+Q2+ as.factor(Q3) + Q21 + as.factor(Q8),
      family = "binomial",
      data = subset(x, country == "PL"),
      x = TRUE)
})
how_many_imputations(mod[['(5)']], cv = 0.05, alpha = 0.05)
mod[['(5)']]  <- mice::pool(mod[['(5)']])

mod[['(6)']] <- lapply(imputed_data, function(x) {
  glm(Q11_7 ~ Q16A + Q16G + Q16H + Q16I + Q16J + as.factor(Q4a) + as.factor(Q7_3) + as.factor(Q7_5) + as.factor(Q14)+ as.factor(Q1)+Q2+ as.factor(Q3) + Q21 + as.factor(Q8),
      family = "binomial",
      data = subset(x, country == "SK"),
      x = TRUE)
})
mod[['(6)']]  <- mice::pool(mod[['(6)']])
how_many_imputations(mod[['(6)']], cv = 0.05, alpha = 0.05)

coef_map_full <- c("Q16H"="Labor Market", "Q16I"="Economy", "Q16A"="Reduce Support" , "Q16G"="Cultural Threat", "Q16J"="Safety", "as.factor(Q7_3)3"="Econ. Situation", "as.factor(Q7_5)5"="Govt. Eval.", "as.factor(Q4a)1"="Urban", "as.factor(Q14)2"="No Inappropriate Beh.", "as.factor(Q1)2"="Female", "Q2"="Age", "as.factor(Q3)2"= "Secondary, no leaving exam", "as.factor(Q3)3"= "Secondary, leaving exam", "as.factor(Q3)4"= "University","Q21"= "Household Financial Status", "as.factor(Q8)2" = "War Responsibility: Ukraine","as.factor(Q8)3" = "War Responsibility: US/NATO", "as.factor(COUNTRY)2"="CZ", "as.factor(COUNTRY)3"="PL", "as.factor(COUNTRY)4"="HU")

modelsummary(model=mod, output="latex", coef_map = coef_map_full, stars = c('*' = .05,'**'=0.01))

rm(list = ls())

#######################################################################################################
#### Study 1: Figure 3: Views on Support for Ukrainian Refugees ####

data <- rio::import("Hlatky_UKR_Ref_Study1_Survey_Data.csv")

data <- data %>%
  mutate(
    Q16A  = case_when(Q16A  == 1 ~ 4,
                      Q16A  == 2 ~ 3,
                      Q16A  == 3 ~ 2,
                      Q16A  == 4 ~ 1,
                      TRUE       ~ Q16A),
    
    Q16C = case_when(Q16C  == 1 ~ 4,
                     Q16C  == 2 ~ 3,
                     Q16C  == 3 ~ 2,
                     Q16C  == 4 ~ 1,
                     TRUE       ~ Q16C),
    
    Q16D = case_when(Q16D  == 1 ~ 4,
                     Q16D  == 2 ~ 3,
                     Q16D  == 3 ~ 2,
                     Q16D  == 4 ~ 1,
                     TRUE       ~ Q16D),
    
    Q16E = case_when(Q16E  == 1 ~ 4,
                     Q16E  == 2 ~ 3,
                     Q16E  == 3 ~ 2,
                     Q16E  == 4 ~ 1,
                     TRUE       ~ Q16E),
    
    Q16H  = case_when(Q16H  == 1 ~ 4,
                      Q16H  == 2 ~ 3,
                      Q16H  == 3 ~ 2,
                      Q16H  == 4 ~ 1,
                      TRUE       ~ Q16H),
    
    Q16J  = case_when(Q16J  == 1 ~ 4,
                      Q16J  == 2 ~ 3,
                      Q16J  == 3 ~ 2,
                      Q16J  == 4 ~ 1,
                      TRUE       ~ Q16J)
  )

data_long <- data %>%
  pivot_longer(cols = c(Q16B, Q16C, Q16D, Q16E), 
               names_to = "item", values_to = "response") %>%
  mutate(
    support = case_when(
      item == "Q16B" ~ ifelse(is.na(response), 0, ifelse(response >= 3, 1, 0)),
      item %in% c("Q16C", "Q16D", "Q16E") ~ ifelse(is.na(response), 0, ifelse(response >= 3, 1, 0)),
      TRUE ~ 0
    )
  )

summary_data <- data_long %>%
  group_by(country, item) %>%
  summarise(
    n_total = n(),                        
    prop_support = mean(support),           
    se = sqrt(prop_support * (1 - prop_support) / n_total),
    .groups = "drop"
  ) %>%
  mutate(
    support_percent = prop_support * 100,
    lower = (prop_support - 1.96 * se) * 100,
    upper = (prop_support + 1.96 * se) * 100
  )

labels <- c(
  Q16B = "Free Healthcare",
  Q16C = "Free Language Courses",
  Q16D = "Free Public Transport",
  Q16E = "Reduced Rent"
)

support_type_plot <- ggplot(summary_data, aes(x = country, y = support_percent, color=country)) +
  geom_point(size = 1.25) +
  scale_color_manual(values = c(
    "CZ" = "#999999",
    "SK" = "#E69F00",
    "PL" = "#56B4E9",  
    "HU" = "#009E73"
  ), guide = "none") +
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0) +
  facet_wrap(~ item, labeller = as_labeller(labels)) +
  labs(
    y = "% Agree",
    x = "Country",
  ) +
  theme_minimal() +
  theme(
    strip.background = element_rect(fill = "#f0f0f0", color = "gray70"),
    panel.border = element_rect(color = "black", fill = NA, linewidth = 0.25),
    strip.text = element_text(size = 10)
  ) + ylim(0,100)

support_type_plot

# ggsave("Fig3_Int_Program_Support.jpeg", plot = support_type_plot, width=6 ,height = 4.5, dpi=600)

#######################################################################################################

#######################################################################################################
#### Study 1: Figure A2: Views on Support for Ukrainian Refugees (Original Scale) ####

data_long <- data %>%
  pivot_longer(cols = c(Q16B, Q16C, Q16D, Q16E), 
               names_to = "item", values_to = "response")

summary_data <- data_long %>%
  filter(!is.na(response)) %>%  # Remove NAs before calculating statistics
  group_by(country, item) %>%
  summarise(
    n_total = n(),                        
    mean_response = mean(response),
    se = sd(response) / sqrt(n_total),
    lower = mean_response - 1.96 * se,
    upper = mean_response + 1.96 * se,
    .groups = "drop"
  )

labels <- c(
  Q16B = "Free Healthcare",
  Q16C = "Free Language Courses",
  Q16D = "Free Public Transport",
  Q16E = "Reduced Rent"
)

support_type_plot_alt <- ggplot(summary_data, aes(x = country, y = mean_response, color = country)) +
  geom_point(size = 1.25) +
  scale_color_manual(values = c(
    "CZ" = "#999999",
    "SK" = "#E69F00",
    "PL" = "#56B4E9",  
    "HU" = "#009E73"
  ), guide = "none") +
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0) +
  facet_wrap(~ item, labeller = as_labeller(labels)) +
  labs(
    y = "Mean Response",
    x = "Country"
  ) +
  theme_minimal() +
  theme(
    strip.background = element_rect(fill = "#f0f0f0", color = "gray70"),
    panel.border = element_rect(color = "black", fill = NA, linewidth = 0.25),
    strip.text = element_text(size = 10)
  )

support_type_plot_alt

# ggsave("FigA2_Int_Program_Support_alt.jpeg", plot = support_type_plot_alt, width=6 ,height = 4.5, dpi=600)

#######################################################################################################

#######################################################################################################
#### Study 1: Figure A3: Social Distance Ukrainian Refugees versus Non-Ukrainian Refugees ####

data_long <- data %>%
  pivot_longer(
    cols = c(Q11_1, Q11_2, Q11_3, Q11_4, Q12_1, Q12_2, Q12_3, Q12_4), 
    names_to = "item", values_to = "response"
  ) %>%
  mutate(
    accept = case_when(
      item %in% c("Q11_1", "Q11_2", "Q11_3", "Q11_4", 
                  "Q12_1", "Q12_2", "Q12_3", "Q12_4") ~ 
        ifelse(is.na(response), 0, ifelse(response == 1, 1, 0)),
      TRUE ~ 0
    )
  )

summary_data <- data_long %>%
  group_by(country, item) %>%
  summarise(
    n_total = n(),                                 
    accept_percent = mean(accept) * 100,          
    se = sqrt((accept_percent / 100) * (1 - accept_percent / 100) / n_total),
    lower = accept_percent - 1.96 * se * 100,
    upper = accept_percent + 1.96 * se * 100,
    .groups = "drop"
  )

summary_data <- summary_data %>%
  mutate(refugee = case_when(
    str_starts(item, "Q11") ~ "Ukrainians",
    str_starts(item, "Q12") ~ "Non-Ukrainians",
    TRUE ~ NA_character_
  ))

labels <- c(
  Q11_1 = "Family",
  Q11_2 = "Close Friends",
  Q11_3 = "Neighbors",
  Q11_4 = "Colleagues",
  Q12_1 = "Family",
  Q12_2 = "Close Friends",
  Q12_3 = "Neighbors",
  Q12_4 = "Colleagues"
)

summary_data$item <- labels[summary_data$item]

summary_data$item <- factor(summary_data$item, 
                            levels = c("Family", "Close Friends", "Neighbors", "Colleagues"))

social_distance_plot <- ggplot(summary_data, aes(x = country, y = accept_percent, shape = refugee, color = country)) +
  geom_point(size = 1.25, position = position_dodge(width = 0.5)) +
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0, position = position_dodge(width = 0.5)) +
  scale_color_manual(values = c(
    "CZ" = "#999999",
    "SK" = "#E69F00",
    "PL" = "#56B4E9",  
    "HU" = "#009E73"
  ), guide = "none") +
  guides(shape = guide_legend(override.aes = list(
    color="#999999",
    size = 1.25
  ), title = NULL,reverse = FALSE)) +
  facet_wrap(~ item) +
  labs(
    y = "% Agree",
    x = "Country",
  ) +
  theme_minimal() +
  theme(
    strip.background = element_rect(fill = "#f0f0f0", color = "gray70"),
    strip.text = element_text(size = 10),
    panel.border = element_rect(color = "black", fill = NA, linewidth = 0.25),
    legend.text = element_text(size = 8),        
    legend.title = element_text(size = 9),       
    legend.key.size = unit(0.6, "lines"),        
    legend.position = "bottom"
  ) + 
  ylim(0,100)

social_distance_plot

# ggsave("FigA3_SocialDist.jpeg", plot = social_distance_plot, width=6, height = 4.5, dpi=600)

rm(list = ls())

#######################################################################################################

#######################################################################################################
#### Study 1: Figure A4: Ukrainian Refugee Admission Preferences and Economic Concerns; Figure A5: Ukrainian Refugee Admission Preferences and Attitudes related to Sociocultural Conservatism ####

eb1 <- rio::import("EB_97_5.dta") 
eb2 <- rio::import("EB_98_2.dta")
eb3 <- rio::import("EB_99_4.dta")
eb4 <- rio::import("EB_100_2.dta")
eb5 <- rio::import("EB_101_3.dta")



# Fix EB1 labels to match EB2
attr(eb1$sd18a, "labels") <- c(
  "Very satisfied"       = 1,
  "Fairly satisfied"     = 2,
  "Not very satisfied"   = 3,
  "Not at all satisfied" = 4,
  "Don't know (SPONT.)"  = 5
)

eb1$year <- "July 2022"
eb2$year <- "Feb. 2023"
eb3$year <- "June 2023"
eb4$year <- "Nov. 2023"
eb5$year <- "May 2024"

eb_all <- bind_rows(
  eb1 %>% transmute(
    uniqid = uniqid, 
    isocntry = isocntry, 
    year = year, 
    household_finances        = qa1_5,
    national_economy          = qa1_2,
    eu_image                  = d78,
    feel_immig_non_eu = NA, 
    accept_ukr_ref = qe2_5
  ),
  eb2 %>% transmute(
    uniqid = uniqid, 
    isocntry = isocntry, 
    year = year, 
    household_finances        = qa1_5,
    national_economy          = qa1_2,
    eu_image                  = d78,
    feel_immig_non_eu = NA, 
    accept_ukr_ref = qe2_5
  ),
  eb3 %>% transmute(
    uniqid = uniqid, 
    isocntry = isocntry, 
    year = year, 
    household_finances        = qa1_5,
    national_economy          = qa1_2,
    eu_image                  = d78,
    feel_immig_non_eu = qb10_2, 
    accept_ukr_ref = qe2_5
  ),
  eb4 %>% transmute(
    uniqid = uniqid, 
    isocntry = isocntry, 
    year = year, 
    household_finances        = qa1_5,
    national_economy          = qa1_2,
    eu_image                  = d78,
    feel_immig_non_eu = qb7_2, 
    accept_ukr_ref = qd2_5
  ),
  eb5 %>% transmute(
    uniqid = uniqid, 
    isocntry = isocntry, 
    year = year, 
    household_finances        = qa1_5,
    national_economy          = qa1_2,
    eu_image                  = d78,
    feel_immig_non_eu = qb8_2, 
    accept_ukr_ref = qd2_5
  )
)


eb_all$household_finances        <- ifelse(eb_all$household_finances %in% c(5), NA, eb_all$household_finances)
eb_all$national_economy          <- ifelse(eb_all$national_economy %in% c(5), NA, eb_all$national_economy)
eb_all$eu_image                  <- ifelse(eb_all$eu_image %in% c(6), NA, eb_all$eu_image)
eb_all$feel_immig_non_eu             <- ifelse(eb_all$feel_immig_non_eu  %in% c(5), NA, eb_all$feel_immig_non_eu)
eb_all <-  eb_all[eb_all$isocntry=="CZ" | eb_all$isocntry=="HU" | eb_all$isocntry=="PL" | eb_all$isocntry=="SK", ]
eb_all$year <- factor(eb_all$year, levels = c("July 2022", "Feb. 2023", "June 2023", "Nov. 2023", "May 2024"))


eb_all <- eb_all %>%
  mutate(accept_ukr_ref = case_when(
    accept_ukr_ref == 1 ~ 4,
    accept_ukr_ref == 2 ~ 3,
    accept_ukr_ref == 3 ~ 2,
    accept_ukr_ref == 4 ~ 1,
    accept_ukr_ref == 5 ~ NA_real_,  # Add this back
    TRUE ~ accept_ukr_ref
  ))

eb_all$feel_immig_non_eu_bin <- ifelse(eb_all$feel_immig_non_eu %in% 1:2, 1, 0)


eb_all <- eb_all %>% 
  mutate(
    
    eu_image_group = case_when(
      eu_image%in% 1:2 ~ "Positive Evaluation",
      eu_image %in% 4:5 ~ "Negative Evaluation", 
      TRUE ~ NA_character_
    ),
    

    household_finances_group = case_when(
      household_finances %in% 1:2 ~ "Positive Evaluation",
      household_finances %in% 3:4 ~ "Negative Evaluation",
      TRUE ~ NA_character_
    ),
    
    national_economy_group = case_when(
      national_economy %in% 1:2 ~ "Positive Evaluation",
      national_economy %in% 3:4 ~ "Negative Evaluation",
      TRUE ~ NA_character_
    ),
    
    feel_immig_non_eu_group = case_when(
      feel_immig_non_eu %in% 1:2 ~ "Positive Evaluation",
      feel_immig_non_eu %in% 3:4 ~ "Negative Evaluation",
      TRUE ~ NA_character_
    )
  )

create_subgroup_data <- function(data, group_var) {
  data %>%
    filter(!is.na(!!sym(group_var))) %>%
    group_by(isocntry, year, !!sym(group_var)) %>%
    summarise(
      n_total = n(),
      n_support = sum(accept_ukr_ref %in% 3:4, na.rm = TRUE),
      pct_support = 100 * n_support / n_total,
      .groups = 'drop'
    ) %>%
    mutate(
      isocntry = case_when(
        isocntry == "CZ" ~ "Czech Rep.",
        isocntry == "HU" ~ "Hungary", 
        isocntry == "PL" ~ "Poland",
        isocntry == "SK" ~ "Slovakia",
        TRUE ~ isocntry
      )
    )
}

hh_finances_data <- create_subgroup_data(eb_all,"household_finances_group")
national_economy_data <- create_subgroup_data(eb_all, "national_economy_group")
eu_image_data <- create_subgroup_data(eb_all, "eu_image_group")
non_eu_immig_data <- create_subgroup_data(eb_all, "feel_immig_non_eu_group")

create_subgroup_plot <- function(data, group_var, title) {
  ggplot(data, aes(x = year, y = pct_support, color = !!sym(group_var), group = !!sym(group_var))) +
    geom_line(size = 0.7) +
    geom_point(size = 1.2) +
    facet_wrap(~ isocntry, scales = "free_y") +
    scale_color_manual(
      values = c("#999999", "#009E73", "#56B4E9")
    ) +
    ylim(20, 100) +
    labs(
      title = paste0(
        "<span style='background-color:#f0f0f0; border:1px solid gray70;
         padding:4px 8px; border-radius:4px;'>", 
        title, "</span>"),
      x = "Survey Date",
      y = "% Agree",
      color = str_to_title(gsub("_", " ", group_var))
    ) +
    theme_minimal() +
    theme(
      plot.title = element_markdown(hjust = 0.5, size = 12, face = "bold"), 
      axis.text.x = element_text(angle = 45, hjust = 1),
      legend.position = "bottom",
      legend.title = element_blank(),
      strip.background = element_rect(fill = "#f0f0f0", color = "gray70"),
      strip.text = element_text(size = 10, face = "bold")
    )
}

household_financial_plot <- create_subgroup_plot(hh_finances_data, "household_finances_group", "Household Financial Situation")
national_economy_plot <- create_subgroup_plot(national_economy_data, "national_economy_group", "National Economy")
eu_image_plot <- create_subgroup_plot(eu_image_data, "eu_image_group", "EU Attitudes")
non_eu_immig_plot <- create_subgroup_plot(non_eu_immig_data, "feel_immig_non_eu_group", "Non-EU Immigration")


household_financial_plot <- household_financial_plot + guides(color = guide_legend(reverse = TRUE)) 
eu_image_plot <- eu_image_plot + guides(color = guide_legend(reverse = TRUE)) 
national_economy_plot <- national_economy_plot  +
  theme(
    axis.title.y = element_blank(),  # remove y-axis title
  ) + guides(color = guide_legend(reverse = TRUE))
non_eu_immig_plot <- non_eu_immig_plot  +
  theme(
    axis.title.y = element_blank(),  # remove y-axis title
  ) + guides(color = guide_legend(reverse = TRUE))


combined_1 <- household_financial_plot + 
  national_economy_plot + 
  plot_layout(nrow = 1, widths = c(1, 1), guides = "collect") & 
  theme(legend.position = "bottom") 

combined_1

# ggsave("FigA4_EB_Econ_Figure.jpeg", combined_1, width = 6.25, height = 5.5, dpi=600)

combined_2 <- eu_image_plot + 
  non_eu_immig_plot + 
  plot_layout(nrow = 1,widths = c(1, 1), guides = "collect") & 
  theme(legend.position = "bottom")

combined_2

# ggsave("FigA5_Sociocultural_Figure.jpeg", combined_2, width = 6.25, height = 5.5, dpi = 600)

#######################################################################################################
